Purpose¶
- Prediction on time series dataset ("CalIt2 Building People Counts Data Set")
- Detect anomalies
- Detect events hours
- Develop yourself with multivariate data
In [1]:
import pandas as pd
import numpy as np
from pandas import read_csv
from pandas import datetime
from matplotlib import pyplot
import holidays
import seaborn as sns
import pandas_profiling
import plotly.offline as py
import plotly.graph_objs as go
import prophet_prediction as prp
py.init_notebook_mode()
pd.set_option('display.max_rows', 500)
pd.set_option('display.max_columns', 500)
pd.set_option('display.width', 1000)
encoding = "utf-8-sig"
In [2]:
#Â Reading Data
df_orig = read_csv('CalIt2.data', header=None, squeeze=True, names = ["flow_id", "date", "time", "count"])
df = df_orig.copy()
df_orig.head(2)
Out[2]:
In [3]:
# extract datetime
date_time = pd.to_datetime(df['date'] + ' ' + df['time'])
# drop date and time to merge as datetime
df.drop('date', axis=1, inplace=True)
df.drop('time', axis=1, inplace=True)
df = pd.concat([df, date_time], axis=1)
df.columns = ['flow_id', 'count', 'datetime']
df.head()
Out[3]:
In [4]:
pandas_profiling.ProfileReport(df)
Out[4]:
Overview
Dataset info
| Number of variables | 3 |
|---|---|
| Number of observations | 10080 |
| Total Missing (%) | 0.0% |
| Total size in memory | 236.3 KiB |
| Average record size in memory | 24.0 B |
Variables types
| Numeric | 1 |
|---|---|
| Categorical | 0 |
| Boolean | 1 |
| Date | 1 |
| Text (Unique) | 0 |
| Rejected | 0 |
| Unsupported | 0 |
Warnings
counthas 5041 / 50.0% zeros Zeros
Variables
count
Numeric
| Distinct count | 55 |
|---|---|
| Unique (%) | 0.5% |
| Missing (%) | 0.0% |
| Missing (n) | 0 |
| Infinite (%) | 0.0% |
| Infinite (n) | 0 |
| Mean | 3.8129 |
|---|---|
| Minimum | 0 |
| Maximum | 62 |
| Zeros (%) | 50.0% |
Quantile statistics
| Minimum | 0 |
|---|---|
| 5-th percentile | 0 |
| Q1 | 0 |
| Median | 0 |
| Q3 | 5 |
| 95-th percentile | 17 |
| Maximum | 62 |
| Range | 62 |
| Interquartile range | 5 |
Descriptive statistics
| Standard deviation | 6.4358 |
|---|---|
| Coef of variation | 1.6879 |
| Kurtosis | 8.3028 |
| Mean | 3.8129 |
| MAD | 4.6912 |
| Skewness | 2.4507 |
| Sum | 38434 |
| Variance | 41.419 |
| Memory size | 78.8 KiB |
| Value | Count | Frequency (%) | |
| 0 | 5041 | 50.0% | |
| 1 | 1042 | 10.3% | |
| 2 | 656 | 6.5% | |
| 3 | 372 | 3.7% | |
| 4 | 255 | 2.5% | |
| 6 | 237 | 2.4% | |
| 8 | 212 | 2.1% | |
| 5 | 212 | 2.1% | |
| 7 | 209 | 2.1% | |
| 9 | 204 | 2.0% | |
| Other values (45) | 1640 | 16.3% |
Minimum 5 values
| Value | Count | Frequency (%) | |
| 0 | 5041 | 50.0% | |
| 1 | 1042 | 10.3% | |
| 2 | 656 | 6.5% | |
| 3 | 372 | 3.7% | |
| 4 | 255 | 2.5% |
Maximum 5 values
| Value | Count | Frequency (%) | |
| 53 | 1 | 0.0% | |
| 54 | 2 | 0.0% | |
| 55 | 1 | 0.0% | |
| 61 | 1 | 0.0% | |
| 62 | 1 | 0.0% |
datetime
Date
| Distinct count | 5040 |
|---|---|
| Unique (%) | 50.0% |
| Missing (%) | 0.0% |
| Missing (n) | 0 |
| Infinite (%) | 0.0% |
| Infinite (n) | 0 |
| Minimum | 2005-07-24 00:00:00 |
|---|---|
| Maximum | 2005-11-05 23:30:00 |
flow_id
Boolean
| Distinct count | 2 |
|---|---|
| Unique (%) | 0.0% |
| Missing (%) | 0.0% |
| Missing (n) | 0 |
| Mean | 8 |
|---|
| 7 | |
|---|---|
| 9 |
| Value | Count | Frequency (%) | |
| 7 | 5040 | 50.0% | |
| 9 | 5040 | 50.0% |
Correlations
Sample
Comments¶
- people counter data between the 24/07/2005 - 05/11-2005 (4months). So there is no seaonality information on data.
- 'counts' data has '0' values %50 frequently
- max value of 'count':62. Most of the values between 0-10. Look at the above profiler code: Variables -> Count -> Histogram
- No missing data
Prepare dataframes¶
'df_in' --> incoming people at main door
'df_out' --> outgoing people at main door
'df_dif' --> incoming - outgoing people at main door
'df_in_out' --> incoming and outgoing people as separete columns
In [5]:
df_in = df.loc[df['flow_id'] == 9]
df_out = df.loc[df['flow_id'] == 7]
df_in.drop('flow_id', inplace=True, axis=1)
df_in.set_index('datetime', inplace=True)
df_out.drop('flow_id', inplace=True, axis=1)
df_out.set_index('datetime', inplace=True)
df_dif = pd.concat([df_in['count'] - df_out['count']], axis=1)
df_in_out = pd.concat([df_in['count'], df_out['count']], axis=1)
df_in_out.columns= ["count_in", "count_out"]
In [6]:
print('df_in')
print(df_in.head(2))
print("--------")
print('df_out')
print(df_out.head(2))
print("--------")
print('df_dif')
print(df_dif.head(2))
print("--------")
print('df_in_out')
print(df_in_out.head(2))
Profiling 'df_in_out'¶
In [7]:
pandas_profiling.ProfileReport(df_in_out)
Out[7]:
Overview
Dataset info
| Number of variables | 3 |
|---|---|
| Number of observations | 5040 |
| Total Missing (%) | 0.0% |
| Total size in memory | 118.2 KiB |
| Average record size in memory | 24.0 B |
Variables types
| Numeric | 2 |
|---|---|
| Categorical | 0 |
| Boolean | 0 |
| Date | 1 |
| Text (Unique) | 0 |
| Rejected | 0 |
| Unsupported | 0 |
Variables
count_in
Numeric
| Distinct count | 48 |
|---|---|
| Unique (%) | 1.0% |
| Missing (%) | 0.0% |
| Missing (n) | 0 |
| Infinite (%) | 0.0% |
| Infinite (n) | 0 |
| Mean | 3.8883 |
|---|---|
| Minimum | 0 |
| Maximum | 62 |
| Zeros (%) | 55.5% |
Quantile statistics
| Minimum | 0 |
|---|---|
| 5-th percentile | 0 |
| Q1 | 0 |
| Median | 0 |
| Q3 | 5.25 |
| 95-th percentile | 18 |
| Maximum | 62 |
| Range | 62 |
| Interquartile range | 5.25 |
Descriptive statistics
| Standard deviation | 6.8346 |
|---|---|
| Coef of variation | 1.7577 |
| Kurtosis | 7.912 |
| Mean | 3.8883 |
| MAD | 5.0044 |
| Skewness | 2.4242 |
| Sum | 19597 |
| Variance | 46.711 |
| Memory size | 39.5 KiB |
| Value | Count | Frequency (%) | |
| 0 | 2799 | 55.5% | |
| 1 | 392 | 7.8% | |
| 2 | 242 | 4.8% | |
| 3 | 156 | 3.1% | |
| 6 | 112 | 2.2% | |
| 8 | 98 | 1.9% | |
| 10 | 98 | 1.9% | |
| 5 | 97 | 1.9% | |
| 11 | 96 | 1.9% | |
| 4 | 94 | 1.9% | |
| Other values (38) | 856 | 17.0% |
Minimum 5 values
| Value | Count | Frequency (%) | |
| 0 | 2799 | 55.5% | |
| 1 | 392 | 7.8% | |
| 2 | 242 | 4.8% | |
| 3 | 156 | 3.1% | |
| 4 | 94 | 1.9% |
Maximum 5 values
| Value | Count | Frequency (%) | |
| 50 | 3 | 0.1% | |
| 53 | 1 | 0.0% | |
| 55 | 1 | 0.0% | |
| 61 | 1 | 0.0% | |
| 62 | 1 | 0.0% |
count_out
Numeric
| Distinct count | 45 |
|---|---|
| Unique (%) | 0.9% |
| Missing (%) | 0.0% |
| Missing (n) | 0 |
| Infinite (%) | 0.0% |
| Infinite (n) | 0 |
| Mean | 3.7375 |
|---|---|
| Minimum | 0 |
| Maximum | 54 |
| Zeros (%) | 44.5% |
Quantile statistics
| Minimum | 0 |
|---|---|
| 5-th percentile | 0 |
| Q1 | 0 |
| Median | 1 |
| Q3 | 5 |
| 95-th percentile | 16 |
| Maximum | 54 |
| Range | 54 |
| Interquartile range | 5 |
Descriptive statistics
| Standard deviation | 6.0104 |
|---|---|
| Coef of variation | 1.6081 |
| Kurtosis | 8.433 |
| Mean | 3.7375 |
| MAD | 4.38 |
| Skewness | 2.4439 |
| Sum | 18837 |
| Variance | 36.125 |
| Memory size | 39.5 KiB |
| Value | Count | Frequency (%) | |
| 0 | 2242 | 44.5% | |
| 1 | 650 | 12.9% | |
| 2 | 414 | 8.2% | |
| 3 | 216 | 4.3% | |
| 4 | 161 | 3.2% | |
| 6 | 125 | 2.5% | |
| 7 | 121 | 2.4% | |
| 5 | 115 | 2.3% | |
| 8 | 114 | 2.3% | |
| 9 | 114 | 2.3% | |
| Other values (35) | 768 | 15.2% |
Minimum 5 values
| Value | Count | Frequency (%) | |
| 0 | 2242 | 44.5% | |
| 1 | 650 | 12.9% | |
| 2 | 414 | 8.2% | |
| 3 | 216 | 4.3% | |
| 4 | 161 | 3.2% |
Maximum 5 values
| Value | Count | Frequency (%) | |
| 43 | 1 | 0.0% | |
| 47 | 1 | 0.0% | |
| 49 | 1 | 0.0% | |
| 52 | 1 | 0.0% | |
| 54 | 2 | 0.0% |
datetime
Date
| Distinct count | 5040 |
|---|---|
| Unique (%) | 100.0% |
| Missing (%) | 0.0% |
| Missing (n) | 0 |
| Infinite (%) | 0.0% |
| Infinite (n) | 0 |
| Minimum | 2005-07-24 00:00:00 |
|---|---|
| Maximum | 2005-11-05 23:30:00 |
Correlations
Sample
Analyze Trend: Showing Yearly, Montly, Daily Trending times for:¶
- 'df_in'
- 'df_out'
- 'df_dif'
'df_in'¶
In [8]:
forecast_in, m_in = prp.prophet_predict(df_in)
mae_in = prp.forecast_mae(forecast_in)
print(mae_in)
In [9]:
# Plot Trend
fig = m_in.plot_components(forecast_in) # Plot
In [10]:
# Plot Forecast
prp.plot_forecast_ts(forecast_in) # Plot
'df_out'¶
In [11]:
forecast_out, m_out = prp.prophet_predict(df_out)
mae_out = prp.forecast_mae(forecast_out)
print(mae_out)
In [12]:
# Plot Trend
fig = m_out.plot_components(forecast_out) # Plot
In [13]:
# Plot Forecast
prp.plot_forecast_ts(forecast_out) # Plot
'df_dif'¶
In [14]:
forecast_dif, m_dif = prp.prophet_predict(df_dif)
mae_dif = prp.forecast_mae(forecast_dif)
print(mae_dif)
In [15]:
# Plot Trend
fig = m_dif.plot_components(forecast_dif) # Plot
In [16]:
# Plot Forecast
prp.plot_forecast_ts(forecast_dif) # Plot
Correlation beetwen incoming and outgoing counter¶
In [17]:
corrmat_in_out = df_in_out.corr().abs()
f, axes = pyplot.subplots(1, 1, figsize=(2, 2))
sns.heatmap(corrmat_in_out, vmax=1, vmin=0, square=True,ax=axes)
corrmat_in_out
Out[17]:
In [18]:
corrmat_in_out
Out[18]:
Regression¶
In [19]:
df_in_regression = df_in.copy()
df_out_regression = df_out.copy()
In [20]:
from datetime import datetime
def make_categorical_hour(hour):
if hour < 4:
return "1"
elif hour < 8:
return "2"
elif hour < 12:
return "3"
elif hour < 16:
return "4"
elif hour < 20:
return "5"
elif hour < 22:
return "6"
elif hour < 25:
return "1"
def add_date_feature(df, df_date):
df['dayofweek'] = df_date.dt.dayofweek.astype(int)
df['hour'] = df_date.dt.hour.astype(int).map(make_categorical_hour)
us_holidays = holidays.UnitedStates(prov=None, state='CA')
df['isholiday'] = [x in us_holidays for x in df_date[date_time.index]]
In [21]:
temp = df_orig.loc[df_orig['flow_id'] == 9]
date_time = pd.to_datetime(temp['date'] + ' ' + temp['time'])
date_time.index = date_time
In [22]:
add_date_feature(df_in_regression, date_time)
add_date_feature(df_out_regression, date_time)
In [23]:
df_in_regression.sample(2)
Out[23]:
In [24]:
df_out_regression.sample(2)
Out[24]:
In [25]:
from sklearn.model_selection import KFold
from sklearn.pipeline import Pipeline
from sklearn.feature_selection import SelectFromModel
from sklearn.metrics import r2_score, explained_variance_score, \
mean_squared_error, mean_absolute_error, \
mean_squared_log_error, median_absolute_error
from xgboost import XGBRegressor
from lightgbm import LGBMRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn import linear_model
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression, Perceptron
from statsmodels.tsa.arima_model import ARIMA
In [26]:
#import pdb
In [27]:
def forecast_cross_validation(df_data, y_name='count'):
y = df_data[y_name]
X = df_data.drop(y_name, axis=1)
kf = KFold(n_splits=7, shuffle=True) # Define the split
Xv = X.values
yv = y.values
preds = []
scores = []
for train_index, test_index in kf.split(X):
X_train, X_test = Xv[train_index], Xv[test_index]
y_train, y_test = yv[train_index], yv[test_index]
pipline_model = Pipeline([
('feature_selection', SelectFromModel(XGBRegressor(), threshold='0.2*mean')),
('classification', XGBRegressor())
])
# Ensemble
xgb = XGBRegressor()
rf = RandomForestRegressor(n_estimators=30)
lgbm = LGBMRegressor()
#Linear Models
reg = linear_model.LinearRegression()
ridge = linear_model.Ridge(alpha=0.5)
lasso = linear_model.Lasso(alpha=0.1)
elasticNet = linear_model.ElasticNet()
omp = linear_model.OrthogonalMatchingPursuit()
sgd = linear_model.SGDRegressor()
polyReg3 = Pipeline([('poly', PolynomialFeatures(degree=2)),
('feature_selection', SelectFromModel(LinearRegression(fit_intercept=False))),
('linear', LinearRegression(fit_intercept=False))])
models = [pipline_model, reg, ridge, lasso, elasticNet, omp, sgd, polyReg3, lgbm, rf, xgb]
pred_df = pd.DataFrame()
score_df = pd.DataFrame()
for model in models:
model.fit(X_train, y_train)
y_pred = model.predict(X_test)
pred_df[type(model).__name__] = y_pred
#pdb.set_trace()
score_df.at[type(model).__name__, "r2"] = r2_score(y_test, y_pred)
score_df.at[type(model).__name__, "rmse"] = np.sqrt(mean_squared_error(y_test, y_pred))
score_df.at[type(model).__name__, "mse"] = mean_squared_error(y_test, y_pred)
score_df.at[type(model).__name__, "mae"] = mean_absolute_error(y_test, y_pred)
score_df.at[type(model).__name__, "exp_var"] = explained_variance_score(y_test, y_pred)
score_df.at[type(model).__name__, "medae_cv"] = median_absolute_error(y_test, y_pred)
scores.append(score_df)
preds.append(pred_df)
n_splits = kf.get_n_splits(X)
#average the splits
scores_average = scores[0]/n_splits
for i in range(1,n_splits):
scores_average =scores_average.add(scores[i]/n_splits, fill_value=0)
scores_average.applymap(lambda x:'{0:.3f}'.format(x))
return scores_average.sort_values(by=['r2'])
In [28]:
forecast_cross_validation(df_in_regression)
Out[28]:
In [29]:
forecast_cross_validation(df_out_regression)
Out[29]:
ARIMA¶
In [32]:
df_in.head(2)
Out[32]:
In [33]:
df_out.head(2)
Out[33]:
In [39]:
df_diff = pd.concat([df_in_out['count_in'] - df_in_out['count_out'] ], axis=1)
df_diff.columns = ["count"]
df_diff.head(3)
Out[39]:
In [52]:
from statsmodels.tsa.arima_model import ARIMA
from tqdm import tqdm
def arima_forecast(df):
X = df.astype('float64').values
size = int(len(X) * 0.66)
train, test = X[0:size], X[size:len(X)]
history = [x for x in train]
predictions = list()
for t in tqdm(range(len(test))):
model = ARIMA(history, order=(5,1,0))
model_fit = model.fit(disp=0)
output = model_fit.forecast()
yhat = output[0]
predictions.append(yhat)
obs = test[t]
history.append(obs)
error = mean_squared_error(test, predictions)
print('Test MSE: %.3f' % error)
pyplot.plot(test, color='blue')
pyplot.plot(predictions, color='yellow')
pyplot.show()
In [53]:
arima_forecast(df_in)
In [54]:
arima_forecast(df_out)
In [55]:
arima_forecast(df_diff)
TODOS¶
In regression:
- add before days (Shifting)
- Add prophet trend informations
- substract incoming people and outgoing people for predictions and look at error
- determine anomalies corresponding to predictions. If error is high there may be anomaly
IN Arıma can be seen the anomalies with graphic
- Try LSTM
- In General:
- Research..
- Entegrate with the 'events' data: build a pipeline to find event start and end hours.
- Combine all of the informations
In [30]:
#from keras.models import Sequential
#from keras.layers import LSTM
#from keras.layers import Dense
## create and fit the LSTM network
#lstm = Sequential()
#lstm.add(LSTM(30))
#lstm.add(Dense(2))
#lstm.compile(loss='mean_squared_error', optimizer='adam')